addpath ('c:\matlab7\fdaM')
addpath ('c:\matlab7\fdaM\examples\weather')

%  Last modified 6 September 2004

%  -----------------------------------------------------------------------
%                     Daily Weather Data
%  -----------------------------------------------------------------------

%  ------------------------  input the data  -----------------------

fid    = fopen('dailtemp.dat','rt');
tempav = fscanf(fid,'%f');
tempav = reshape(tempav, [365,35]);

fid    = fopen('dailprec.dat','rt');
precav = fscanf(fid,'%f');
precav = reshape(precav, [365,35]);

%  set up the times of observation at noon

daytime = (1:365)'-0.5;

%  day values roughly in weeks

weeks = linspace(0,365,53)';   

%  define 8-character names for stations

place = [ ...
    'Arvida     '; ...
    'Bagottville'; ...
    'Calgary    '; ...
    'Charlottown'; ...
    'Churchill  '; ...
    'Dawson     '; ...
    'Edmonton   '; ...
    'Fredericton'; ...
    'Halifax    '; ...
    'Inuvik     '; ...
    'Iqaluit    '; ...
    'Kamloops   '; ...
    'London     '; ...
    'Montreal   '; ...
    'Ottawa     '; ...
    'Pr. Albert '; ...
    'Pr. George '; ...
    'Pr. Rupert '; ...
    'Quebec     '; ...
    'Regina     '; ...
    'Resolute   '; ...
    'Scheffervll'; ...
    'Sherbrooke '; ...
    'St. Johns  '; ...
    'Sydney     '; ...
    'The Pas    '; ...
    'Thunder Bay'; ...
    'Toronto    '; ...
    'Uranium Cty'; ...
    'Vancouver  '; ...
    'Victoria   '; ...
    'Whitehorse '; ...
    'Winnipeg   '; ...
    'Yarmouth   '; ...
    'Yellowknife'];

%  set up indices that order the stations from east to west to north

geogindex = [24,  9, 25, 34,  4,  8, 22,  1,  2, 19, 23, 14, 15, 28, 13, ...
             27, 33, 26,  5, 20, 16, 29,  7,  3, 12, 30, 31, 17, 18, 32, ...
              6, 35, 11, 10, 21];
         
%  set up the harmonic acceleration operator

Lbasis  = create_constant_basis([0,365]);  %  create a constant basis
Lcoef   = [0,(2*pi/365)^2,0];    %  set up three coefficients
wfd     = fd(Lcoef,Lbasis);      % define an FD object for weight functions
wfdcell = fd2cell(wfd);          % convert the FD object to a cell object
harmaccelLfd = Lfd(3, wfdcell);  %  define the operator object

% set up a saturated basis capable of interpolating the data

nbasis      = 365;  
daybasis365 = create_fourier_basis([0,365], nbasis);
wtvec       = ones(365,1);

lambda = 1e7;  %  minimum GCV estimate, corresponding to 8.5 df
fdParobj = fdPar(daybasis365, harmaccelLfd, lambda);

%  set up some quantities

daybasis365  = getbasis(getfd(fdParobj));
Zmat         = eval_basis(daytime, daybasis365);
harmaccelLfd = getLfd(fdParobj);
Rmat         = eval_penalty(daybasis365, harmaccelLfd);

%  put the stations in geographical order

tempav = tempav(:,geogindex);
precav = precav(:,geogindex);
place  = place(geogindex,:);

%  use interactive menu program to try out different plots

station =   1;  %  number of St. John's
lambda  = 100;  %  a smoothing level that corresponds to 56 df
[H_f1, H_f2, Hc_OK, ...
      Hc_Next, Hc_Last, Hc_This, Hc_Enter, Hc_Quit,  ...
      Hc_Weather, Hc_Velocity, Hc_Accel, Hc_HarmAcc, Hc_Residual, ...
      Hc_Temp, Hc_Prec, Hc_CaseNo, Hc_Lambda] = ...
      dailysetup(station, lambda, daytime, ...
                 tempav, precav, ...
                 fdParobj, Zmat, Rmat, place);
pause;
delete(H_f1);
delete(H_f2);



